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Abstract 

Numerical  solutions  have  been  obtained  for  steady  natural  convec- 
tion in  a square  cavity.  The  numerical  method  used  was  developed  for 
unsteady,  incompressible,  viscous  fluid  flow.  The  similarity 
parameters  were  chosen  to  match  those  of  an  international  comparison 
exercise.  Results  are  presented  and  compared  with  those  obtained  by 
other  researchers  using  different  methods. 

Key  Words:  Cavity  flow;  fluid  dynamics;  natural  convection;  numerical 

methods . 

Introducti  on 

The  idea  of  carrying  out  a comparison  exercise  for  numerical 
schemes  based  on  the  problem  of  steady  natural  convection  in  a square 
cavity  with  vertical  sides  which  are  maintained  at  different 
temperatures  was  proposed  in  1979  [1].  The  problem  is  of  interest  due 
to  the  range  of  detailed  flow  patterns  found  in  this  simple  geometry. 
The  applicable  equations  here  are  the  Navier-Stokes  and  continuity 
equations  of  fluid  flow.  The  Navier-Stokes  equations  include  a body 
force  term,  which  requires  the  solution  of  an  additional  equation  for 
temperature.  These  equations  are  all  simplified  by  the  Boussinesq 
approximation.  There  are  no  singularities  in  the  boundary  conditions. 
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as  occur  in  the  driven  cavity  problem.  The  results  of  this  study  could 
be  used  as  a model  problem  solution  for  checking  computer  codes 
developed  for  practical  elliptic  fluid  flow  computations. 

Two  recent  papers  summarize  the  contributions  received  from 
participating  researchers  [2]  and  a bench  mark  solution  [3]  with  which 
their  results  were  compared.  This  report  describes  a solution  obtained 
at  NBS  using  a numerical  scheme  developed  for  unsteady,  incompressible, 
viscous  fluid  flow  problems. 

In  the  following  section,  the  comparison  problem  and  similarity 
parameters  will  be  presented.  The  solution  procedure  will  then  be 
described.  Results  and  comparisons  with  the  bench  mark  solution  and 
those  of  other  contributors  are  also  included. 


The  Comparison  Problem 

The  basic  problem  considered  here  is  the  steady  two-dimensional 
flow  of  a Boussinesq  fluid  in  an  upright  square  cavity.  The  details 
are  shown  in  figure  1.  Each  side  of  the  cavity  has  length  D,  and  both 
velocity  components  are  zero  on  the  walls.  The  vertical  sides  are  at 
temperatures  T-|  and  J2>  with  the  normal  temperature  gradient  being 
zero  on  the  top  and  bottom  walls.  The  fluid  has  thermal  diffusivity  < 
and  kinematic  viscosity  v,  with  Prandtl  number  Pr  =^.  The  Rayleigh 

3 < 

number  for  this  flow  is  defined  as  Ra  = where  3 is  the  coeffi- 

<v 

cient  of  expansion  of  the  fluid,  g is  gravitational  acceleration,  and 
At  = Ti  ■ T£-  In  this  study,  all  lengths  are  nondimensional ized 

K 

with  respect  to  D;  all  velocities  with  respect  to  jj ; time  with  respect 
D2 

to  — ; and  p,  the  ratio  of  pressure  to  reference  density,  with  respect 

*2  T'  - T2 

to  ^75-.  The  nondimensional  temperature  is  T = f f—  , where 

D2  * 1 " ‘ 2 
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T'  is  the  dimensional  value.  Therefore,  0 <_  x <_  1 , 0 <_  z _<  1 , T = 1 at 
x = 0 and  T = 0 at  x = 1 . 

The  velocity  and  temperature  fields  are  to  be  computed  with 
Prandtl  number  0.71  for  Rayleigh  numbers  of  103,  104,  10^,  and 
10^.  The  comparison  also  includes  Nusselt  numbers,  or  nondimensional 
values  of  heat  flux  in  the  horizontal  direction. 


Numerical  Modeling 

The  two-dimensional  Navier-Stokes , continuity,  and  temperature 
equations  for  a viscous  Boussinesq  fluid  are 
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where  u and  w are  velocity  components  in  the  x-  and  z-  directions, 
respecti  vely. 

The  basic  numerical  scheme  employed  here  uses  quadratic  upwind 
differencing  for  convection  and  an  explicit  Leith-type  of  temporal 
differencing  [4].  A fast  direct  method  is  used  to  solve  the  Poisson 
equation  for  pressure  at  each  time  step.  The  solution  is  accomplished 
on  a staggered  mesh  in  which  pressures  and  temperatures  are  defined  at 
cell  centers  and  normal  velocities  at  cell  faces. 
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The  boundary  conditions  on  the  four  cavity  walls  are  that  normal 

and  tangential  velocities  are  zero.  The  temperatures  on  the  vertical 

walls  are  specified.  Any  values  for  grid  points  just  outside  the  walls 

are  obtained  by  quadratic  extrapolation  using  the  wall  value  and  two 

points  just  inside.  The  extrapolated  values  are  used  in  computing 

3T 

derivatives  at  the  wall.  At  the  top  and  bottom  walls,^-  = 0,  so  the 
points  just  outside  those  walls  have  the  same  temperatures  as  those 
just  inside  the  walls. 

The  only  mesh  employed  in  this  study  was  a 50  x 50  uniform  grid 
with  Ax  = Az  = 0.02.  The  initial  conditions  for  the  computations  were 
either  u = 0,  w = 0,  T = 1-x  for  Ra  = 103,  or  the  results  of  a 
previous  computation  at  a lower  Rayleigh  number.  The  time  step,  At, 
was  set  at  a maximum  of  7.5  • 10~3  so  that  the  diffusion 
coefficient  in  the  temperature  equation  was  less  than  0.5  and  the 
Courant  number  less  than  one.  Computation  times  on  the  MBS  UNIVAC 
1100/82  required  to  obtain  steady-state  solutions  ranged  between  5 and 
7 hours.  At  steady-state,  the  averages  of  the  absolute  values  of  the 

time-derivatives  of  the  velocities  [- — ^ and  3 — , where  n 

indicates  time  level]  throughout  the  cavity  were  0(10“3)  for  the 
lower  three  Rayleigh  numbers  and  0(10""*)  for  Ra  = 10°.  The 
maximum  change  in  u or  w during  the  final  time  step  was  0(10“^)  for 
Ra  = 10°  and  0(10~3)  for  the  lowest  Rayleigh  numbers.  The 
average  of  the  absolute  values  of  the  time-derivatives  of  temperature 
was  0(10"3)  or  less.  The  summation  of  the  absolute  values  of  the 
mass  flux  residuals  was  0(10-3)  for  Ra  = 103  and  increased  by  an 
order  of  magnitude  for  each  tenfold  increase  in  Ra. 
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Results  and  Discussion 


Computations  have  been  carried  out  for  the  four  Rayleigh  numbers 
in  the  comparison  exercise.  The  results  are  shown  in  table  I.  The 
quantities  listed  are  those  requested  for  comparison  purposes: 

umax»  z - the  maximum  u-component  of  velocity  along  the  line 
x = 0.5  and  its  location; 

wmax»  x ” the  maximum  w-component  of  velocity  along  the  line 
z = 0.5  and  its  location; 

Nu0  - the  average  Mussel t number  on  the  wall  of  the  cavity  at 
x = 0; 

Mumax,  z " the  maximum  value  of  the  local  Mussel t number  on 
the  wall  at  x = 0 and  its  location; 

Mumi-n,  z - the  minimum  value  of  the  local  Mussel t number  on 
the  wall  at  x = 0 and  its  location; 

Mu  - the  average  Mussel t number  throughout  the  cavity; 

Mui/2  " the  average  Mussel t number  along  the  line  x = 0.5. 

The  values  for  umax,  z,  wmax,  and  x are  grid  point 

3T 

values  (no  interpolation).  The  Musselt  number,  Mu(x,z)  = uT  was 

calculated  at  each  of  the  mesh  points.  The  quantities  Munax  and 
Mumi-n  and  their  locations  are  values  of  Nu(0,  z)  at  the  z grid 
points.  The  average  Nusselt  number  along  a line  x = constant,  such  as 

fl 

Mu  or  Mu, /?,  was  obtained  by  integration:  Mu  = J Mu(x,z)  dz. 

o M c x - o 

The  integrals  were  evaluated  by  the  rectangular  rule.  Average  values 
obtained  at  each  value  of  x were  integrated  by  the  trapezoidal  rule  to 
find  Mu  = j Muxdx. 

Contour  plots  of  temperature  T,  velocity  components  u and  w,  and 
vorticity  are  presented  in  figures  2-5  for  each  of  the  four  Rayleigh 
numbers.  These  illustrate  the  diagonal  symmetry  of  the  problem: 
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T(  x , z)  = 1 - T( 1 -x,  1-z),  u(x,z)  = -u(l-x,  1-z),  w(x,z)  = -w(l-x,  1-z). 
Also,  the  development  of  the  thermal  boundary  layer  as  Ra  increases  is 
clearly  seen.  The  same  qualitative  features  appear  in  the  bench  mark 
solution  (figures  4-7  in  reference  3). 

Figure  6 presents  velocity  vector  plots  for  the  four  cases,  where 
the  line  segments  indicate  the  local  flow  direction.  They  are 
analogous  to  the  streamline  plots  in  figure  3 of  reference  3. 

Secondary  recirculating  rolls  appear  for  Ra  j>  10^  as  observed 
prevdoti-s-ly  [1,5]. 

The  quantitative  results  presented  here  in  table  I can  be  compared 
with  those  of  the  bench  mark  solution  given  in  table  I of  reference  2. 
The  percentage  differences  (with  respect  to  the  bench  mark  solution) 
are  shown  here  in  table  II.  The  field  locations  where  the  parameters 
of  table  II  occur  have  been  omitted  from  the  comparisons  as  was  done  in 
reference  2.  Also,  stream  function  was  not  computed  in  the  present 
work.  The  percentage  differences  are  seen  to  be  small,  within  about  1% 
for  the  lower  three  Rayleigh  numbers. 

It  is  interesting  to  compare  these  results  with  those  in  tables 

V 1 1 - X of  reference  2 which  contain  the  contributed  results.  For 

purposes  of  comparison,  the  absolute  values  of  the  percentage  errors  in 

Mu,  Nu  , Mu  . , u , and  w were  averaged  for  each  solution  at 
* max’  min’  max  ’ max  3 

each  Rayleigh  number.  Mu  is  either  Mu  or  Mu-]/2»  if  available, 
and  has  been  compared  with  Mu-j/2  1,1  the  bench  ^ark  solution.  The 
errors  at  each  Rayleigh  number  were  then  averaged  to  give  an  overall 
average  error  for  each  contributed  solution. 

The  average  errors  in  the  present  work  ranged  from  0.1%  at  Ra  = 

103  to  2.1%  at  Ra  = 10^.  The  error  averaged  over  the  four  Rayleigh 
numbers  is  0.3%.  Mot  all  contributors  computed  all  quantities  or 
submitted  solutions  for  all  Rayleigh  numbers.  Of  the  26  solutions  from 


a 


reference  2 that  had  all  quantities,  four  had  a smaller  average  error 
than  the  present  work,  placing  these  results  in  the  top  20%  of  the 
solutions. 

Conclusion 

Solutions  have  been  presented  for  laminar  natural  convection  in  a 
square  cavity.  The  computer  code  was  adapted  from  one  used  for 
unsteady,  incompressible,  viscous  fluid  flow  problems.  The  results 
obtained  here  compare  well  with  the  bench  mark  solution  given  in 
reference  3.  The  present  results  also  compare  well  with  the 
contributed  solutions  in  reference  2.  Although  the  errors  in  the 
present  solution  increased  with  increasing  Rayleigh  number,  it  should 
be  noted  that  all  four  solutions  were  obtained  using  the  same  grid. 

Mesh  refinement  would  have  been  preferable  as  Ra  increased  but  was  not 
possible  due  to  limitations  on  computer  resources. 
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Table  I.  Results 


Ra 


ioi 

IQ4 

105 

ioi 

umax 

3.646 

16.126 

34.62 

64.86 

z 

0.81 

0.83 

0.85 

0.85 

w 

max 

3.689 

19.507 

67.98 

211.86 

X 

0.17 

0.11 

0.07 

0.03 

Nu0 

1.118 

2.244 

4.526 

8.902 

%ax 

1.507 

3.538 

7.812 

18.809 

Z 

0.09 

0.15 

0.07 

0.03 

Nu  . 
min 

0.691 

0.586 

0.732 

0.997 

z 

0.99 

0.99 

0.99 

0.99 

ni f 

1.118 

2.244 

4.523 

8.893 

Mu1/2 

1.118 

2.244 

4.524 

8.898 

Table 

II.  Percentage  Differences 

Ra 

103 

IQ4 

5 

10° 

Nu0 

0.1 

0.3 

0.4 

1.0 

NlT 

0.0 

0.0 

0.1 

1.1 

Nu1/2 

0.0 

0.0 

0.1 

i.i 

Numax 

0.1 

0.3 

1.2 

4.9 

Nu  . 
mi  n 

-0.1 

0.0 

0.4 

0.8 

umax 

-0.1 

-0.3 

-0.3 

0.4 

w 

max 

-0.2 

-0.6 

-0.9 

-3.4 
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